#ifndef MATRIX_HPP
#define MATRIX_HPP
//=========================================================
#include "vector.hpp"
//=========================================================
//  CLASS MATRIX NxN
//=========================================================
class matrixes{
  size_t dim;
  double **mtrx;
public:
  matrixes(){
    this->dim = 0;
    this->mtrx = NULL;
  };
  //-------------------------------------------------------
  matrixes(size_t size){
    if(size != 0){
    this->dim = size;
    this->mtrx = new double*[this->dim];
    for(size_t i=0; i<this->dim; i++)
      this->mtrx[i] = new double[this->dim];
    for(size_t i=0; i<this->dim; i++)
      for(size_t j=0; j<this->dim; j++)
        this->mtrx[i][j] = 0;
    }
  };
  //-------------------------------------------------------
  matrixes(const matrixes &m){
    this->dim = m.dim;
    if(m.mtrx == NULL){
      this->mtrx = NULL;
      return;
    }
    this->mtrx = new double*[this->dim];
    for(size_t i=0; i<this->dim; i++)
      this->mtrx[i] = new double[this->dim];
    for(size_t i=0; i<this->dim; i++)
      for(size_t j=0; j<this->dim; j++)
        this->mtrx[i][j] = m.mtrx[i][j];
  };
  //-------------------------------------------------------
  ~matrixes(){
    if(this->mtrx != NULL){
      for(size_t i=0; i<this->dim; i++)
        delete [] this->mtrx[i];
      delete []this->mtrx;
    }
  };
  //-------------------------------------------------------
  size_t ret_dim() const{
    return this->dim;
  };
  //-------------------------------------------------------
  //-------------------------------------------------------
  double* operator[](const size_t indx){
    return this->mtrx[indx];
  };
  //-------------------------------------------------------
  const double* operator[](const size_t indx) const{
    return this->mtrx[indx];
  };
  //-------------------------------------------------------
  //-------------------------------------------------------
  //matrix = matrix
  const matrixes& operator=(const matrixes &rhs){
    if(this != &rhs){
      if(this->dim != rhs.dim){
        for(size_t i=0; i<this->dim; i++)
          delete [] this->mtrx[i];
        delete [] this->mtrx;
        this->dim = rhs.dim;
        this->mtrx = new double*[this->dim];
        for(size_t i=0; i<this->dim; i++)
          this->mtrx[i] = new double[this->dim];
      }
      for(size_t i=0; i<this->dim; i++)
        for(size_t j=0; j<this->dim; j++)
          this->mtrx[i][j] = rhs.mtrx[i][j];
    }
    return *this;
  };
  //-------------------------------------------------------
  //matrix * matrix
  matrixes operator*(const matrixes &m){
    matrixes r(this->dim);
    if(this->dim == m.ret_dim()){
      for(size_t i=0; i<this->dim; i++)
        for(size_t j=0; j<this->dim; j++)
          for(size_t k=0; k<this->dim; k++)
            r[i][j] = r[i][j] + this->mtrx[i][k]*m[k][j];
    }
    return r;
  };
  //-------------------------------------------------------
  //matrix * vector
  const matrixes& operator*(const vector &rhs){
    matrixes m(this->dim);
    if(this->dim == rhs.retdim()){
      for(size_t i=0; i<this->dim; i++)
        for(size_t j=0; j<this->dim; j++)
          for(size_t k=0; k<this->dim; k++)
            m[i][j] = m[i][j] + this->mtrx[i][k]*rhs[j];
    }
    *this = m;
    return *this;
  };
  //-------------------------------------------------------
  //vector * matrix
  friend vector operator*(const vector &lhs, const matrixes &m){
    size_t size = m.dim;
    vector c(size);
    if(lhs.retdim() == m.dim){
      for(size_t i=0; i<size; i++)
        for(size_t j=0; j<size; j++)
          c[i] = c[i]+lhs[i]*m.mtrx[j][i];
    }
    return c;
  };
};
//---------------------------------------------------------
//---------------------------------------------------------
//print matrix
void print_matrix(const matrixes &p,FILE *fi){
  size_t dim = p.ret_dim();
  if(dim != 0){
    for(size_t i=0; i<dim; i++){
      for (size_t j=0; j<dim; j++)
        fprintf(fi,"   %.4f",p[i][j]);
      fprintf(fi,"\n");
    }
  }else{
    fprintf(fi,"dimension of the matrix = 0\n");
  }
}
//-------------------------------------------------------
//inverse matrix
matrixes inv_matrix(const matrixes &m){
  double c;
  size_t size = m.ret_dim();
  matrixes e(size), m_obr(size);
  for(size_t i=0; i<size; i++)
    e[i][i] = 1;
  m_obr = m;
  for(size_t k=0; k<size-1; ++k){
    for(size_t i=k+1; i<size; ++i){ 
      c = m_obr[i][k];
      for(size_t j=k; j<size; ++j){
        for(size_t l=0; l<size; ++l){
          e[i][l] = e[i][l]-(m_obr[i][k]/m_obr[k][k])*e[k][l];
        }
        m_obr[i][j] = m_obr[i][j]-(c/m_obr[k][k])*m_obr[k][j];
      }
    }
  }
  for(size_t k=1; k<size; ++k){
    for(size_t i=k-1; i>0; --i){
      c = m_obr[i][k];
      for(size_t j=k; j<size; ++j){
        for(size_t l=0; l<size; ++l){
          e[i][l] = e[i][l]-(m_obr[i][k]/m_obr[k][k])*e[k][l];
        }
        m_obr[i][j] = m_obr[i][j]-(c/m_obr[k][k])*m_obr[k][j];
      }
    }
  }
  for(size_t i=0; i<size; i++){
    c = m_obr[i][i];
    for(size_t j=0; j<size; j++){
      e[i][j] = e[i][j]/c;
      m_obr[i][j] = m_obr[i][j]/c;
    }
  }
  return e;
}
//=========================================================
#endif